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ABSTRACT 


‘This study considers the accuracy of the finite difference method 
in the solution of linear elasticity problems that involve either a stress 
discontinuity or a stress singularity. Solutions to three elasticity 
problems are discussed in detail: a serai-infinite plane subjected to a 

uniform load over a portion of its boundary; a bimetallic plate under 
uniform tensile stress; and a long, midplane symmetric, fiber-reinforced 
laminate subjected to uniform axial strain. 

Finite difference solutions to the three problems are compared 
with finite element solutions to corresponding problems. For the first 
problem a comparison with the exact solution is also made. 

The finite difference formulations for the three problems are based 
on second order finite difference formulas that provide for variable 
spacings in two perpendicular direct j-cns , Forward and backward 
difference formulaj are used near boundaries where their use eliminates 
the need for fictitious grid points. Moreover, forward and backward 
finite difference formulas are used to enforce continuity of interlaminar 
stress components for the third problem. 

The study shows that the finite difference method employed 
in this investigation provides solutions to the three elasticity problems 
considered that are as accurate ns the corresponding finite element 
solutions. Furthermore, the finite difference method appears to give 
a solution for the laminate problem that characterizes the stress 
distributions near an interface corner in a more realistic manner than 
the finite element method. 
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X. INTRODUCTION 


A serious failure mechanism for laminated composite materials is edge 
delamination. Various numerical methods have teen used in attempts to 
calculate the interlaminar stress components that accompany delamination in 
a finite-width [+ 45 ]g angle-ply laminate under uniform axial strain 
[1,2, 3, 4]. These efforts have resulted in serious discrepancies in reported 
behavior for the interlaminar normal stress distribution near an interface 
corner [4J. For example, a finite-difference procedure [1] and a 
perturbation procedure [2] predict tensile Interlaminar normal stress near 
an interface corner, while finite element methods [3,4] predict compressive 
normal stress in this region. Furthermore, some uncertainty exists 
regarding the character of the In-plane, Interlaminar normal and shearing 
stress distributions near an interface corner that are predicted by finite 
element methods [3,4]. 

The primary purpose of this investigation is to determine if the finite 
difference method is capable of providing accurate predictions for Che 
interlaminar stress components near an interface corner and, hence nea a 
stress singularity. 

A second purpose of this investigation is to determine if predictions, 
by finite element methods, for in-plane, interlaminar stress components near 
an interface corner accurately represent laminate behavior; and, if these 
predictions are spurious, to cast light on the origin of the weakness in the 
finite element method that results in the spurious behavior. 

In this investigation the finite difference method has been used to 
obtain numerical solutions for three different problems that involve 
a point where a stress component becomes discontinuous or singular. These 
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problems are: (a) uniform pressure on part of a semi-infinite plane 

(Figure la) , (b) a bimetallic plate under uniform axial tension (Figure 5a) , 
and (c) a finite-width [+ 45] angle-ply laminate under unifona axial 

S 

strain (Figure 8a) , Solutions to each of these problems via finite element 
methods are reported in reference [AX-' 

The finite difference procedure used in this investigation provides for 
variable grid spacings in two perpendicular directions. Consequently, 
computational efficiency is effected by taking closely spaced grid lines in 
regions where the stress components are expected to vary rapidly, and a 
coarser grid in regions where the stress components do not vary rapidly. 

The coefficient matrix corresponding to the system equations is 
unsymmetricalj therefore, it is necessary to store the entire band of the 
coefficient matrix, Moreover, an equation solver capable of handling un- 
symmetrical systems of algebraic equations must be available. Nevertheless, 
variable grid capability leads to more efficient computations than finite 
diffi rence procedures that use uniform spacing because substantially fewer 
grid lines are needed to realize an accuracy comparable to the accuracy 
associated with a specific uniform grid. 


II. DISTRIBUTED LOAD ON A SEMI-INFINITE PLANE 

Figure la depicts a semi-infinite plane that is sub j "^cted to a uniform 

pressure on part of the edge y=0, The exact solution for this problem is given 

in reference [4] and indicates that = + p/'n' when the points (+ a,0) 

are approached along the lines x = +a. Consequently, 0^^ ^ at these points. 

It is of interest in this investigation to obtain a numerical solution for 

this problem based on the finite difference method, and to compare the finite 

difference, finite element, and exact solutions for the stress distributions 

(a , 0 , a ) along the lines x - f a, 
x' y xy — 
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For discretization purposes it is assumed that the stress component are 
nearly zero for x >_ + 10a, and that vertical displacements are cissentially zero 
at a depth y ^ 10a. Moreover, for computational efficiency use is made of 
symmetry with respect to the line x » 0. 

Boundary Value Problem . The field equations associated with the distributed 
load problem are listed below. 


Stress-Strain Relations 


)RIGINAL 


1-v 

E 


(e “t* ve ) 
2 y' 


1-v 

as (J 


2 ^''y + 


xy yx 2(l-fv) xy 
Strain-Displacement Relations: 


E - u 
X ,x 

£ “V 

y »y 


p S3 p =* U + V 

V “,y ^ '^.x 


Equilibrium Equations (Plane Stress) : 


(u 


vv 

,xx 


»xy 

(v 

+ 

vu _ ' 

>yy 


,xy 


.yy 


,xx 


,xy' 

'U ) 
jxy' 


Boundary Conditions: 


u(0,y) = V „(0,y) = 0 
a^(10a,y) = o^y(10a,y) = 0 


^ s y s 


Oy^(x,0) - 0, 0 J< X _< loa; cjy(x,o) = ( -p/2 


-p 0 _< X < a 
X = a 
0 a < X < 10a 


a^^(x,10a) = v(x,10a) = 0 


0 < X < iGa 


( 1 ) 


( 2 ) 


(3) 


( 4 ) 
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Numerical Results. The finite-difference grid used to analyze this problem 


is shown in Figure lb. VtM’t.ical grid lines are more closely spaced on either 
side of the line x * a, while the horizontal grid lines are more closely 
spaced near the line y ■ 0. Numerical values of grid spacings for the 
X and y directions are listed in Table I. The numerical results to be 
discussed are based on these spacings which correspond to 2,146 degrees of 
freedom. 

la Figures 2 and 3 the open circles and dashed lines represent numerical 
solutions obtained via the finite difference and finite element methods, 
respectively, and the solid lines represent the exact solution, for the 
stress components and along the line x » a. The finite dif- 

ference and finite element solutions for o^(a,y) and Oy(a,y) exhibit excellent 
agreement everywhere. The finite difference and finite element solutions 

for 0 (a,y) show excellent agreement with the exact solution except near 

xy 

the point (a,0) where the finite difference solution appears to provide a 

somewhat better agreement - except for the first two nodes of the finite 

difference grid. The finite difference solution for a (a,y) is "drawn” 

xy 

to zero by the enforced zero shearing stress at the boundary, while the 
finite element solution is "drawn" down but not to a zero value at the 
boundary. 

It appears that requiring the stress tensor to be symmetrical at the 
point (a,0) affects the finite difference solution for the shearing stress 
CTj^y.(a,y) only in a small region that is confined to the first two finite 
difference grid points. This region can be made as small as desired, con- 
tingent on numerical limitations. 




The finite difference solution for the stresses o , and cr shown 

X y xy 

in Figure 2 is based on a boundary value 0^ ■ ~ p /2 at the point (a, 0 ) 5 

that is, on an average of the boundary load intensity to the left and right 

of the point (a, 0 ) . A finite difference solution using 0y(a,O) «* - p diffet^ij: 

from this solution only in a small region near the point (a,0) as shown by 

the open circles and solid curves in Figures 4 a, 4 b, and 4 c. From these 

figures it is seen t1*at 0j^y(a,y) is essentially the same for either 

cf„(a,0) * - p/2 or a (a,0) * - p, while the solutions for 0 and 0 are 
j y y X 

affected dramatically in the vicinity of y » 0 . Otherwise, the finite 
difference solutions using 0^(a,O) » - p /2 or 0y(a,O) » - p are essentially 
identical. 

Since the stress tensor is unsymmetrical at the point (a, 0 ) it was of 

interest to determine if a more accurate representation of the behavior of 

0j^(a,y) could be obtained near the point (a,0) by discarding the symmetry 

re,l«itlvi'i: 0 « 0 at this point and replacing it with a finite moment 

xy yx 

equation that would require 0 (a, 0 ) « 0 , but o„^(a, 0 ) 5^ 0^„(a,O). In 

yA yx xy 

addition to the finite moment equation;, a finite force equilibrium was intro- 
duced. The stress distributions for 0^ and 0^ are essentially identical to 
the finite-difference solution for 0y(a,O) = - p/2 and 0yj^(a,O) = o^y(a,0) *0 
everywhere (solid lines in Figures 4 b and 4 c). The shearing stress distribu- 
tion differs only in the neighborhood of the point (a, 0 ). The shearing stress 
0,5.^(a,y) for the case a (a,0) ^ a (a,0) is indicated by the dashed line in 
Figure 4 a. 
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III. BIMETALLIC PUTE UNDER UNIFORM TENSION 


Figure 5a depicts a bimetallic plate under uniform », ensile stress along 
the edges y * + with stress-free boundaries at the edges x « 0, 8a. A 
numerical solution for the stress components along the bond line, based on 
the finite element method, is given in reference [4] for a rigid bottom plate. 

It is of interest in this investigation to obtain a numerical solution 
for the stress components along the bond line using the finite-difference 
method, and to compare this solu'ion with the finite element solution obtained 
in reference [4j. It is of particular interest to observe whether the finite 
difference method is capable of predicting the behavior of the shearing stress 
component near the intersections of the bond line and the free edges . 

Boundary Value Problem . The plane strain field equations for the bimetallic 
plate with a rigid bottom plate are obtained from the plane stress field 
equations given by Equations (l)-(3) by replacing E and v in these equations 
by E = E/(l-v ) and v* = v/(l-v) and affixing the boundary conditions 


(7j^(0>y) = a^y(o,y) = 0 

u(4a,y) = V (4a, y) = 0 
» 


u(x,0) = v(x,0) - 0 

o^^(x,8a)=0, cTy(x,8a)=p 


0 5 y £ 8a ' 


0 £ X £ 4a ) 


(5) 


These boundary conditions make use of symmetry with respect to the plate 
centerline x = 4a. 
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Numerical Results . Figure 5b shows the finite-difference grid used to analyze 
the bimetallic plate problem. Since the bond line stress components are 
expected to change rapidly near the singular point 0, the finite difference 
grid lines are more closely spaced in the region near point 0. Numerical 
values of spaclngs for the x and y directions are listed in Table II. 

Numerical results presented in this section are based on these grid spacings 
which correspond to 2,340 degrees of freedom. 

Corner points of a rectangular finite difference grid are usually 
troublesome because a decision must be made as to which of two possible sets 
of boundary conditions to employ there. In the present investigation it was 
physically appealing to require the displacement components (u,v) at the 
corners of the bond line to be specified (as zero) , since the two plates do 
not separate there. Moreover, at the left corner of the loaded edge (0,8a) 
boundary conditions associated with the stress free edge were imposed, 
while at the right corner of the loaded edge (4a, 8a) the conditions 
cr^ = p and u=0 were imposed. Boundary conditions and grid points to which 
they apply are shown in Figure 5b. 

The open circles and dashed curves in Figure 6a represent the finite 
difference and finite element predictions for the shearing stress 
distribution along the bond line, respectively. This figure indicates that 
the finite difference method has no more trouble predicting the stress 
distribution along the bond line than the finite element method. Indeed the 
two numerical solutions are essentially the same. 
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Figure 6b shows the shearing stress and normal stress distributions along 


the free edge x “ 0 based on the finite difference method. Both a and a 

xy X 

are zero at every finite difference grid point at which these stresses were 

required to be zero. They were nonzero only at the corner point of the bond 

line. It is noted that values of o and o at this corner point are calculated 

xy X 

from the stress-strain relations and represent limiting values of internal 
stresses as the corner point is approached along the bond line. They are not 
necessarily the boundary values on the edge x = 0 at y = 0. This observation 
again suggests that the stress tensor is unsyrametrical at a stress singularity. 

Figure 7 shows a comparison of the bond line shearing stress distribution 
for two different finite difference grids. The solid curve with open circles 
represents the finite difference prediction based on the grid spacings shown 
in Table II. This curve is an exploded view of the behavior or the shearing 
stress o^^CKjO) near the point 0 that is exiiibited in Figure 6a. The dashed 
curve with open squares represents the finite difference prediction based on 
the grid spacings shown in Table III. This finite difference grid maintains 
the same number of rows and columns as the grid of Table II, but the grid lines 
parallel to both the x and y directions are redistributed so that they are more 
dense near point 0. 
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IV. FOUR-PLY LAMINATE UNDER UNIFORM AXIAL STRAIN 

Figure 8a depicts a long, xnidplane symmetric laminate of width 2b. The 

laminate consists of four plies, each of thickness h, and is loaded by a 

uniform axial strain e . Various numerical methods have been used by 

o ■ 

different investigators [1,2, 3, 4,] to predict the distributions of normal 
and shearing stresses between adjacent lamina. Of particular importance is 
the reliability of a particular numerical method to provide a reasonably 
accurate assessment of the behavior of tlje interlaminar stress components 
near the intersection of an interface with a free edge. This point of 
Intersection is referred to as the interface corner [4] and is shown in 
Figure 8a. 

Computations based on the finite element method have yielded stress 
distributions that appear to be reasonable for all interlaminar stress 
components except very near the ii..terface corner. At the interface corner 
the predicted distributions for the inplane, interlaminar stress components 
(o^jOy, and cr^y) tend to digress from a logical extrapolation of the 
stress distributions predicted for interior points along the interface. It 
is of Interest in this investigation to determine whether this digressive 
behavior exhibited by the finite element method represents actual laminate 
behavior or, if the predictions are spurious, to illuminate the origin of 
the weakness in the finite element method that results in this spurious 
behavior. 

A second objective of this investigation is to assess the viability of 
the finite difference method as an effective numerical method in the 
computation of interlaminar stress distributions, particularly near an 
interface corner. 
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ORIGINAL PAGE IS 
OF POOR QUALITY 


It is customary when dealing with this problem to make use of geometric 
and material symmetries, thereby making it necessary to consider only the 
part of the laminate that lies in the first quadrant of the yz plane. This 
part of the laminate is emphasized by the cross-hatched area in Figure 8a. 
The heavy dot in this figure is at the interface corner. 

Boundary Value Problem . The field equations [1] associated with the 
four-ply, [+ laminate are listed below. 

Stress-Strain Relations: 


=' V „ + C-o W + C,, U 

X 11 0 12 ,y 13 ,z — 16 ,y 

y 12 0 11 ,y 13 ,z — 16 ,y 


z 13' 0 33 ,z — 36 ,y 


= C,,(W , + V ) 


zy "44 '‘",y 


( 6 ) 


a = C,/ U 
zx 44 ,z 


V=i"l6^"0 + ^y> ±S6 ^z-^C^6 ^y 


Equilibrium Equations : 


± 

+ 


°66 'J.yy + Ss ± ^26 ^yy ± ^36 «,y. = ° 

<=26 “,yy + =22 ^yy + ^4 + «=23 + =44>".yz = ° 

“,yz + «=23 + «445\yz + <=44 ",yy + <=33 ”,zz = ° 


( 7 ) 


Displacement distribution: 

u(x,y,z) = Eq X + U(y,z) '' 
v(x,y,z) = V(y,z) 
w(x,y,z) = W(y,z) ^ 


( 8 ) 
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Boundary Conditions ; 


ORIGINAU PAGE 5S 
OF POOR QUALITY 


U(0,z) « V(0,z) » W y(0,z) » 0 
U(0,0) « V(0,0) = W(0,0) » 0 


0 < z £ 2h 
z = 0 


(9a) 

(9b) 


Oyj(b,z) = cr^^^(b,z) - ‘^y^g,(b,z) ■= 0 

OyJb,z) = ay^^(b,2) = ‘’y^^(b,z> = 0 
U ^(y,0) = V (y,0) = W(y,0) = 0 

, , <5 

o,u<y>2W = “zxu<y-2h) = = 0 


0 < z < h 


(9c) 


h < z < 2h 


(9d) 


0 £ y £ b [ 


0 £ y £ b 

/ 


(9e) 

(9f) 


(9g) 


In Equations (6) and (7) the upper sign (plus sign) is associated with the 
upper ply (+ 45 ply) and the lower sign (minus sign) is associated with the lower 
ply (- 45 ply) . Equations (8) are fundamental assumptions regarding the dis- 
tribution of the displacement components u,v and w and are given in reference [1] . 

Equations (9a) are the conditions associated with laminate symmetry with 
respect to the z axis, and Equations (9b) are required to exclude rigid body 
motions. Equations (9c) and (9d) require that the edge at y = b be stress-free, 
except at the interface corner. Equations (9e) result from symmetry conditions 
with respect to the y axis, and Equations (9f) require that the edge at z = 2h 
be stress-free. Finally, Equations (9g) require that the interlaminar stress 
components be continuous across the interface. 
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It is particularly important to observe that along the stress-free boundary 
(y “ b) the formulas that express the stress components 0 ^ and in terms 
of displacements (Equations 6) are different for the upper and lower p] ies . 
Consequently, the boundary conditions that should be applied at the point (b,h) 
are not immediately obvious. This observation is a possible clue as regards the 
behavior of the finite element method near the interface corner. 

Numerical Results . The following strategy was used to formulate a finite dif- 
ference model of the four-ply laminate. 

Initially each ply is considered to occupy a separate, independent region. 
Separate, independent finite difference grids are assigned to the regions 
occupied by the two plies. Subsequently, the finite difference module cor- 
responding to the equilibrium equations that are associated with a particular 
ply is applied to each grid point that does not lie on the boundary of the 
region occupied by that ply. The two regions are connected appropriately by 
requiring the displacements (U,V,W) and the interlaminar stress components 

(o , a , and o ) be continuous across the boundary common to the two plies, 

2 2y 2X 

This approach leads logically to the required boundary condition at the 
interface corner, That is, the interlaminar stress components should be 
required to be continuous across the interface at the interface corner. Thus, 
the need to formulate a boundary condition at the interface corner that 
accounts for the boundary stresses associated with the +45 and -45 plies in 
an equitable manner is avoided. 

It will be observed later that the preceding strategy results ita predic- 
tions for stress distributions that agree well with the finite element predic- 
tions away from the interface corner, and also behave in a much more logical 
manner near the interface corner. Furthermore, the affect that prescribing 
boundary stresses at the interface corner, (instead of interlaminar stress 


12 - 


continuity) has on the stress distributions will be demonstrated. These observa- 
tions provide a clue as to the puzzling behavior of the finite element method 
near the interface corner. 

The finite difference grid used to analyze this problem is shown in Figure 9. 
Since the displacement components (U,V,W) are required to be continuous across 
the interface, the grids associated with the two plys are shown connected in 
Figure 9. Vertical grid lines are more closely spaced near the interface corner 
where the stress components are expected to change rapidly, and the horizontal 
grid lines are more closed spaced about the interface. Numerical values of grid 
spacings for the y and z directions are listed in Table IV. The numerical 
results to be discussed are based on these spacings which correspond to 1989 
degrees of freedom. 

Boundary conditions and the grid points to which they are applied, for what 
is referred to here as the principal finite difference solution, are shown in 
Figure 9. 

Finite difference solutions for three other sets of boundary conditions 
at the interface corner are also discussed. These solutions require either 


Q =(j =0 =0 or a n ^ a a „ = 0 or (a) =(a ) 

yu yxu yzu yxH yzZ y ave yx' ave 


(a ) = 0 

yz ave 


at the interface corner. All other boundary conditions remain the same for 


each of the four cases. Here = (a + a»)/2 at the interface corner. 

ave u <0 

Because the interlaminar stress continuity requirements at the interface 


corner must be relaxed when any of the described sets of conditions are employed, 
discontinuities in the interlaminar stress components (a , o , and a ) are 

Z Z Y 

expected at the interface corner. 


Case I; Principal Finite Difference Solution . As was stated previously the 
boundary conditions and the grid points to which they are applied are depicted 
in Figure 9. Especially important is that continuity of the interlaminar stress 
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components (a^, cr^_, and is required along the interface, including the 

interface corner. 

Figures 10 and 11 compare the finite element and the finite difference 
predictions for stress distributions along the interface between the +45 and 
-45 plies. The open circles connected by a solid curve represent the finite 
difference predictions and the open squares connected by dashed curves re- 
present the finite element solution. 

It is convenient in discussing the behavior of the stresses along the 

interface to segregate them into two groups; the in-plane components a , a , 

X y 

and a , and the interlaminar components o , a , and a . The first group of 
xy’ z’ ZX' zy 

stresses must be identified with a particular ply, even at the interface, be- 
cause they are calculate! from stress-strain relations that are different for 
each ply. The second group of stresses act between the plies and are truly 
interlaminar stresses. They are equal in magnitude owing to the interlaminar 
stress continuity requirements ^ 

Consider first the interlaminar stress components. Figure 10 indicates 

that the finite difference and finite element predictions for the normal stress 

a are in excellent agreement. Most importantly both predict a large compressive 
z 

stress at the interface corner and a small tensile region just interior to the 
interface corner. The lower most curve in Figure 11 depicts the finite dif- 
ference distribution for a along the interface. The finite element prediction 

zx 

essentially coincides with the finite difference prediction except at the inter- 
face corner and is not shown in the figure. The two numerical methods do, how- 
ever, predict similar behavior at the interface corner; that is, the existence 
of a stress singularity for the component o . 

The finite difference method predicts a = 0 along the interface, 

zy 

including the interface corner. This agrees with the finite element prediction 
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except near the interface corner where the finite element method predicts a 
.udden increase in 

Now consider the distributions of the in-plane stress components (a„, o ) 

A y Ay 

along the interface. Figure 11 compares these distributions with corresponding 
distributions predicted by the finite element method. 

Figure 11 shows that the finite element and finite difference predictions 
for each of the in-plane stress components are in excellent agreement except near 
the interface corner. The finite difference method suggests that the in-plane 
stress components become singular at the interface corner, while the finite 
element method predicts a sudden attenuation in the stress components at the 
interface corner. 

Figure 12 shows a comparison of the finite difference and finite element 
predictions for the distribution of a along the free edge. Excellent agreement 
is again observed. 

Figures 13, 14 and 15 show the variation in free edge stress components 

a (b,z), a (b,z), and a Ah,z) as reported in reference [4]. The present finite 
y yz yx 

difference predictions show that a^(b,z), Oy^(b,z) are identically zero every- 
where along the free edge except at the interface corner, and that ay^(b,z) is 
zero everywhere, including the interface corner. It should be noted that 

a (b,z), a (b,z), and a (b,z) are required by the finite difference method to 
y yx yz 

be zero at all grid points along y = b except the grid point that coincides xvrith 
the interface comer. Thus, the values for these latter stress components that 
are shown in the figures are calculated from the stress-strain relations and may 
not represent boundary values. That is the stress tensor may not be symmetrical 
at the interface corner. 
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It is known that at a singular point all stress components are either zero 
or are singular with the same power. The finite difference solution presented 
here appears to satisfy this criterion. 

Case II, a * a » a » 0 at The Interface Corner, The Case II finite 

■ ■■Mll. fi .i w I III! 

difference solution differs from the principal finite difference solution 
only in the boundary conditions applied at the grid point that coincides with 
the interface corner. Accordingly, Interlaminar stress continuity at the inter- 
face corner is replaced by specifying that a (b,h) «= u (b,h) ® <7 _ (b,h) ■ 0, 

yu yxu yzu 

That is, the stress components on the free edge that are associated with the 
upper ply are prescribed to be zero at the interface corner . 

Figure 16 shows the finite difference predictions for the distributions 
of the interlaminar stress components along the interface. This figure shows 
that discontinuities in the normal stress and the shearing stress 
occur at the interface corner. Otherwise, continuity of the interlaminar group 
is maintained along the interface. 

Figure 17 shows the finite difference predictions for the distributions 
of the in-plane group of stress components along the interface for the upper 
and lower plies. An important observation to make from the curves in Figure 17 
is that each stress component associated with the lower ply 
behaves as if a stress singularity existed at the Interface corner, while each 
stress component associated with the upper ply ^y^ *^xyu^ shows a sudden 

and drastic digression from what appears to be a distribution that is trying to 
follow the corresponding distribution of the lower ply. 

The stress distributions associated with the in-plane stresses of the 
upper ply exhibit behaviors near the interface corner similar to those exhibited 
by the finite element method. 
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Figure 18 shows the distribution of the boundary stresses along the free 

edge y » b« The finite difference method requires Oy^(b>z) ■ Oyj^^(b, 2 ) » 

a (b,z) •« 0 at all grid points in the upper ply, including the interface 
y 2U 

corner. It requires that a^^(,btz) « " 0 at all grid 

points in the lower ply, except at the interface corner. Numerical values for 
cJy^(b,h), Oyj^j^(b,h), and calculated from the appropriate stress- 

strain relations are shown on the figure. 


Case III, cr ^ « ^yzi * ^ Interface Corner. The Case III finite 

difference solution differs from the principal finite difference solution only 
in the boundary conditions applied at the grid point that coincides with the inter- 
face corner. Accordingly, interlaminar stress continuity at the interface 
corner is replaced by the conditions a^^(b,h) « 

That is. the stress components on the free edge that are associated with the 
lower ply are prescribed to be zero at the interface corner. 


Figure 19 shows the finite difference distributions for the interlaminar 


stress components. The stress component a^y(y,h) is sero everywhere along the 

interface. This figure shows that discontinuities in the normal stress a and 

z 

the shearing stress o occur at the interface in exactly the same manner as 

Z A 

in Case II. 


Figure 20 shows the finite difference predictions for the distributions of 
the in-plane group of stress components along the interface. An important ob- 
servation to make from the curves in Figure 20 is that each stress component 
associated with the upper ply behaves as if a stress singularity 

existed at the interface corner, while each stress component associated with the 
lower ply cr^^, cr^yjj^) shows a sudden and drastic digression from what 

appears to be a distribution that is trying to follow the corresponding distribu- 
tion of the upper ply. 
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Figure 21 shows the distributions of the boundary stresses along the free edge 

y-b. The finite difference method requires -0 at 

all grid points in the lower ply, including the Interface corner. It requires 

c7„,,(b ,z) * o._,,(b,z) •» 0 (b,z)*0 at all grid points in the upper ply, except 

yu yxu yzix 

at the Interface corner. Numerical values for 0 . (b,h), 0 „ (b,h), and 0 ,. ..(b,h) 

yu ' yxu yzu 

calculated from the appropriate stress-strain relations are shown on the figure. 

Case IV. Average Stress Boundary Condition. The Case IV finite difference 
solution differs from the principal finite difference solution only in the boundary 
conditions imposed at the grid point that coincides with the interface corner. 
Accordingly, interlaminar stress continuity at the interface corner is replaced 
by the conditions “ ^^yx^ave " ^^yz^ave * ® interface corner. 

Here (0)„„_ denotes the average of corresponding stress components for the upper 
and lower plies, 

Figure 22 shows the finite difference distributions for the interlaminar 
stress components along the interface. The finite difference predictions for the 
interlaminar stress components (0 , 0 , and 0 ) exhibited in this figure are 

Z ZX 

in excellent agreement in every respect with the finite element predictions for 
the corresponding stress components. 

Figure 23 shows the finite difference predictions for the distributions of 

the in-plane group of stress components along the interface for the upper and 

lower plies. An important observation to make from the curves in Figure 23 is that 

the stress components associated with both the upper ply (0„,,, i,» 

xui xyu 

the lower ply cr^yil’ *^yP show the curious digressive behavior at the 

interface corner that characterizes the finite element predictions for these 
stress components. 

The foregoing observations provide a clue to the reason the finite 
element, method behaves in the curious manner described in reference [4] near the 
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interface corner. The poaaible reason for the curious behavior is discussed 
in the conclusions section of this report. 

Figure 24 shows the distributions of the boundary stresses along the 

free edge y " b. The finite difference method requires a„„(b,z) •» cr„„„(b, 2 ) » 

yu yxu 

•" 0 at all points in the upper ply, except at the interface corner, 
end Cfy^(b,z) ■ o^y^ih^z) » CJy 2 j^(b,z) - 0 at all points in the lower ply, 
except at the interface corner. Numerical values for cfyy(b,h), n^j^^jCbjh) 
and cTyj 2 u(b,h), Oyj^(b,h), Oyj^j^(b,h), and Cy^^ihth) calculated from appropriate 
stress-strain relations are shown on the figure. Again, it is noted that 
these numerical values are not necessarily boundary values because the 
stress tensor may not be symmetrical at the interface corner. 

An interesting feature of the present case is that while continuity of 
the interlaminar stress components was not enforced directly at the interface 
corner, continuity of these stresses occurred there never thelet.s. 


V. CONCLUSIONS 


Distributed Load Problem . Both the finite element and finite difference 
methods predict shearing stress distributions along the llnex«a that are 
in excellent agreement with each other and with the exact solution, except 
near the point (a,0). The finite difference solution behaves somewhat 
better than the finite element solution in this region, deviating from the 
exact solution only to satisfy the imposed boundary condition Oy^(a,0) » 0, 
Since the exact solution shows that the stress tensor is not symmetrical 
at the point (a,0) , this deviation from the exact solution cannot be 
attributed to an inherent weakness of the finite difference method. The 
trouble arises because of the need to specify a limiting value for o^y(a,0) 
at a point where a a . 

y A 

Bimetallic Plate Probl em. The finite difference and the finite element 

methods predict essentially the same shearing stress distribution along the 

bond line. It is interesting to observe that the displacement components 

are prescribed along the bond line, including the sii.gular point. 

Consequently, no finite difference boundary conditions are prescribed for 

a and a at the grid point that coincides with the singular point. 

X xy 

Numerical values for o and o at the singular point are calculated from 

X 

the stress -strain relations and are not necessarily boundary values since 
the stress tensor may not be symmetrical there. 

Four-Ply Laminate . Based on the numerical evidence presented in the 
principal finite difference solution for the four-ply laminate under 
uniform axial strain it appears that there is no inherent weakness in the 
finite difference method that prevents it from providing accurate 
predictions for the distributions of the interlaminar stress group and for 
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the in-plane stvess group along the interface, except very near the interface 
corner. And, even near the interface corner, the finite difference method 
provides solutions that behave in a way that is characteristic of the behavior 
of stresses near a singular point. 

The four finite difference solutions for the stress distributions associated 
with the interlaminar stress group (a , a) and the in-plane stress group 

Oy, differ only near the interface corner. This is not unexpected, 

since the finite difference models for the four solutions differ only in the 
boundary conditions imposed at a single boundary point-boundary conditions 
that are, moreover, very similar. 

Of the four finite difference solutions the stress distributions predicted by 
the Case I model (stress continuity along the interface) behave near the inter- 
face corner as one expects them to behave near a stress singularity. Therefore, 
it is felt that the Case I predictions should be the definitive solution. 

The Case II and III models exhibit behaviors t^,ut are similar near the 
interface corner. That is, discontinuities appear in the interlaminar stress 
components a and a , and corresponding in-plane stress components for the 
+ 45 ply and the -45 ply show diverger t behaviors near the interface corners. 
Specifically, for the Case II model the in-plane stress components for the +45 
ply are "drawn" down to satisfy the boundary conditions imposed at the inter- 
face corner, while the in-plane stress components for the -45 ply appear to grow 
unboundedly. Just the reverse is true for the Case III finite difference model. 

Stress distributions predicted by the Case IV finite difference model agree 
extremely well, in all respects, with the finite element predictions of the 
same distributions. Consequently, one is led to examine the boundary conditions 
imposed by the two models at the interface corner to explain the curious behavior 
exhibited there by the finite element model. 
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Consider two finite elements, located at the interface corner, which 
Phare the interface corner as a common node. Let one side of each element lie 
on the free edge. The finite element procedure replaces a distributed load on 
an element boundary with concentrated forces acting at the nodes of the element. 
Statical equivalency between the distributed boundary load and the concentrated 
nodal forces is maintained by requiring the virtual work of the nodal forces on 
the corresponding nodal displacements be equal to the virtual work of the actual 
boundary load distribution on the displacements along the boundary to which 
the load is applied. Therefore, the finite element node coincident with the 
interface corner receives "average” contributions from the finite elements on 
either side of the interface. In a finite element solution a stress-free boundary 
condition translates into a nodal force-free boundary condition. Therefore, 
setting the nodal force at the interface corner equal to zero is, in some 
sense, an averaging procedure similar to the boundary conditions used in the 
Case IV model. It is this "averaging" process that apparently eliminates the 
detail that an accurate solution near the interface corner requires. 
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Table IV. Finite-difference grid spacings for the 
four-ply laminate. The grid contained 17 rows and 
39 columns (1989 degrees of freedom) 


Number 
of spaces 

First 

16 

Next 

2 

Next 

4 

Next 

16 

y-spacing 

1.00 

0.50 

__J 

0.25 



0.125 




Number 
of spaces 

First 

4 

Next 

8 

Next 

4 


z-spacing 

0.20 

0.05 

0.20 

. 1 


- 25 - 





























0.5 0.4 0.3 0.2 0.1 0.0 -0.1 -0.2 -0.3 -0.4 



o 




S’ 


c 

ca 


rt 


o 

4-1 

M 

b 

o 

•ri 

o 

w 

■u 

O (U 
cd w 
« 

0) r-i 

a 

'b 

c <u 

to u 


O 4-^ 
B ff 


H 
rH I 
iU -H 

, s 

(1) <u 
4J W 
•H 

b " 
•H T3 
44 0) 
T3 
44 CO 
O O 


cO 
■H 

a 44 
E U 
O to 
u p, 


CO 

(U 

b 

tiO 

*H 

P4 










oniGINAl B 
OF POOR QUALITY 


(b) 


1.0 


2.0 


—— — L 

C 

c 

0 

cr 

o 

1 -T 1 T 

8l 

PO 

*• C 

\ Moment equation result ““ 

o 0 y(asO)« -p 

Q follows -p/2 solution. 

a^(a,o)= -p/2 

\ 

[ 1 1 1 Aj 1 I 1 


- 1.0 

(c) 


0.0 


1.0 


- 0.8 


- 0.6 


y(a,y)/p 
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Figure 4 (continued) Comparison between <?^(a,y), cTy(a,y), 

and a (ajy) for various boundary conditions at the point 
xy 

(a,o) and the special moment equation procedure. 
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(b) Variable spacing finite difference grid. 


Figure 5. Problem involving a stress singularity. 
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0 12 3 4 


X 

a 

(a) Shear stress distribution along the bond line* (y=0) . 



(b) Stress distributions 



/p anda^/p) along the line x=0. 


Figure 6. Cormarison of finite difference and finite element solutions 
for the bimetallic plate problem. 
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26 X 45 finite difference grid spacing 
given in Table II. 


,-n.- 26 X 45 finite difference grid spacing 
given in Table III. 



Comparison of shearing stress distribution along the 

bond line for two finite difference grids. 


(a) 


Figure 


Four-ply laminate 


(b) 3-D stress components 


8. Laminate configuration, loading, and stresses. 
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a "0 "0 *0 

z zx zy 



u -V -v^o 

,z ,•! 


a , a , and a are required to be continuous 
z zx zy 

across the interface except at left interface comer. 


Figure 9. Variable spacing finite-difference grid for the four-ply 
laminate . 











°y(y,h) 

MPa 



(b-y)/h 


Figure 11- Stress discribution along the inter-^ace 
when interlaminar stress continuity is required at 
every point along the interface, including the 
interface corner. 
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I^itcomb, et,al (Fine mesh) 



Figure 13. Distribution of a through the thickness at the free edge* in the +A5 ply 







Figure 14. Distribution of a through the thickness at the free edge in the +45 degree ply. 





through the thlckenss at the free edge in the +45° ply. 
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Figure 16. Stress distributions for the interlaminar 
stress group for a ^ a == a = 0 at the interface 

yLl yzLi 


corner. 








Figure 18. Stress distributions along the boundary v=b for 0 =a =0 = 

° ® ' yu yxu yzu 

the Interface comer. 



a(y,h) 

MPa 



Figure 19. Stress distributions for the interlaminar 
stress group for Oy£ = = <^yzZ - 0 at 

the interface corner. 
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(b-y)/h 


Figure 20. Stress distributions for the in-plane stress 
group for Oyjj = = 0 at the inter 

face corner. 
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Figure 22. Stress distributions for the interlaminar 

stress group for (a ) ave = (a ) ave = (a ) ave = 0 at 

y yx yz 

the interface corner. 
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(b-y)/h 


Figure 23. Stress distributions for the in-plane stress 

group for (a ) ave ® (a ) ave = (a ) ave = 0 at the 
y yx yz 

interface corner. 
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APPENDIX A. FINITE-DIFFERENCE CONSIDERATIONS 


Because the stress components associated with each of the three 
problems considered in this investigation change rapidly only in the 
vicinity of a stress discontinuity or a stress singularity, it is 
computationally effective to use a finite-difference grid with variable 
spacing. Finite-difference formulas for first and second order derivatives 
are easily derived using appropriate Taylor series. Forward, backward, 
and central difference formulas for first order derivatives* as well as 
central difference formulas for second order derivatives, are listed here 
for convenience. Figure A1 depicts essential parameters that appear in 
these equations. 

The finite-difference approximations for first order derivatives are: 


^3^i+2 

(Forward) 

(Al) 

^^l^i-2 ^2^i-l ^3^i 

(Backward) 

(A2) 

^i“ ^l^i-1 ^3^i+l 

(Central) 

(A3) 



A prime is used to denote a first order derivative, l^/hen h and k in these 
formulas are replaced with t and m, respectively, first order derivatives 
for a perpendicular direction are obtained. 

The central finite-difference approximations for second order 
derivatives are : 
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where 


f "1 

j-X 

^i-1' j 
^l-l..i+l 
j-X 
3 

iti._ 

^i+X, j-X 
^i+X, j 
^i+X, j+X 


^ 

^X “ h(h+k) 

2 

^2 “ ~ hk 
2 

^3 “ k(h+k) 


^4“ 

z2 

^5 -&1 

= ^ 

^6 ~m(£-hn) 


km 

^7 ” h,^(h+k) (£+m) 

k (£-m) 

^8 hiia(h+k) 

-k£ 

^9 hm(h+k) (£+m) 


= ni(h-k) 
^XO “hk£(£+m) 


~hm 

^3 ”k^(h+k)(£-hn) 


- (h-k) (^-m) 
XX " hk£m 


-h(£-m) 

^X4 '^"F^Ck+h) 


« = -(h-k)£ 

^X2 hkm(Bmi) 


^ hi 

^15 km (h+k) (Z+m) 


(A5) 


(A6) 


j”l J J-M 



X.r 1 

X 

X “ 1 


(a) Two-dimensionaX centraX 
differences . 


^ X+z 



(b) Forward differences 




i> ■ i ■ ^ 



(c) Backward differences 








-k-Ji 


Figure AX. Finite difference grid 
notations . 


(d) CentraX differences 
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As an example of the derivation of these formulas consider the 

central difference, finite difference formulas. The Taylor series 

expansion of a function f(x) in the neighborhood of a point x^ is 

f(x + Ax) “ f(x ) *f f(x ) Ax + ijjfCx )(Ax)^ + 0{(Ax)^}. 

0 0 0 o 

Now using Figure Al-d write 


and 



< . M 2 

= f, + f^k + hf, 
11 1 

/ . " 2 

= f - f^h + hf,h 
1 i i 


(A7) 


(A8) 


Simultaneous solution of these algebraic equations for f^ and f^ yield the 
appropriate formulas given in A3 and A5. Replacing h and k by £ and m, 
respectively, gives the corresponding finite difference formulas for the 
perpendicular direction. Only the formula for the mixed derivative 
remains to be determined. 

To derive a mixed second order derivative write 

^i,j ^lx^i-l,j '^2x^i,j '^Sx^i+ljj (A9) 

for the line j of Figure Al-a. By differentiation 

^ / • i • 


f. . = Ct f. - . + C„ f . , -f C- f, . 

i,j lx i-l,j 2x i ,3 3x i+l,j. 


(AlO) 


Now 


■i-l,j = ''ly^i-l,j-l "^2y^i-l,j %^i-l,j+l 


(All) 


and similarly for f , and f . , Substituting these formulas into 

1 , j 1 ' 1 j J • 

Equation (AlO) leads to a formula that expresses the mixed derivative 

f,, as a linear combination of the nine nodal points that surround point 
ij 

i,j. The literal subscripts in the coefficients c^^ and c^^ indicate 
these coefficients are to be evaluated using spacings in the x or y 
directions as the subscript dictates. 
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APPENDIX B 


Special Material Element for Semi-Infinite Plane Problem, 


The exact elasticity solution for the partially loaded, semi-infinite 

plane reveals that the shearing stress tends to + p/ir as the 

points (+ a,o) are approached along the lines x - + a(see Figure la) . 

However, the shearing stress, a (+ a,o) applied to the external 

yx — , 

boundary is zero. Consequently, it appears that the stress tensor is not 
symmetrical at these points. 

The equilibrium equations on which the finite difference solution is 

based incorporate symmetry of the stress tensor at all interior points. 

Moreover, no special expressions exist that define the relationship between 

shearing stress components at points where the stress tensor is not 

symmetrical. Therefore, the classical finite difference procedure can not 

be expected to detect such an anomaly. 

Part of this investigation involves examining the effectiveness of 

introducing a moment equation for a finite element of material near point 

(a,o) . The purpose of the moment equation is to establish a relationship 

among the shearing stresses of an unsymmetrical stress tensor. 

Supposedly, the moment equation associated with a finite element of 

material near the boundary point (a,o) replaces a finite difference boundary 

condition involving the shearing stress a (a,o) Problems arise, however, 

yx . 

because once the finite element of material has been introduced it is 
inappropriate to ignore force equilibrium of the element. Thus, three 
equilibrium equations are obtained to replace two finite difference boundary 
conditions at point (a,o). Since no new independent variables have been 
introduced, the resulting system of equations is over specified. 
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Even if it were possible to introduce an additional, appropriate, 
independent stress variable at point (a,o) the system of equations, when 
expressed in terms of displacements, would remain over specified unless 
unsymmetrical stress-strain relations are introduced. This latter concept 
represents a considerable complication and is ignored in the analyses made 
in this investigation. 

Not withstanding the analytical difficulties enumerated in the 
preceding paragraphs, two separate finite elements of material near point 
(a,o) were considered. These elements are shown in Figures B1 and B2. 

The element shown in Figure B1 is based on first order Lagrange 
interpolations of the stresses along its edges. Moment equilibrium for 
the material element yields the equation 


m 


-m 


(m-y)a (a,y)dy- 




Jo 


m 


(m-y)a^(a-£, y)dy- 


Jo 


(x-a+£.)a^(x,ra) dx 


= 0. (Bl) 


Using the first order Lagrange interpolations for stresses along the edge 
of the element leads to the algebraic equation 


m 


(2o , - 2cr „ +a „ -a .) 
x4 x2 x3 xl 


nib / N 

i ^^xy4 xy3^ 






0 . 


(B2) 


Here o a . , a . are components of stress at the finite difference nodes 
XI yi xyi 

that coincide with the corners of the element. 
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Equation B2 is expressed in terms of displacements (u,v) by means 
of the s 3 Tnmetrical stress-strain relations (This is not a strictly 
legitimate procedure since, if the stress tensor is assumed to be 
unsymmetrical, the strain tensor must also be unsymraetrical) . Equation B2 
expressed in terms of displacements is not recorded here because nothing 
of significant interest can be extracted from it. 

Using only the moment equation for the material element and a finite 

difference normal stress boundary condition, encouraging, but not exact, 

agreement with the exact solution was observed for all three stress 

components (a a a ) , 

X, y, xy' 

Because two of the three equilibrium equations associated with the 
material element were ignored it is quite possible that whatever agreement 
that exists, is, probably, gratuitous. 

The second material element is shown in Figure B2 and is based on 
second order Lagrange interpolations of stresses along its edges. Thus, 
approximations of stresses along the edges of the element are of the 
same order of approximation as the finite difference formulas used to solve 
the equilibrium equations. It was felt that this material element 
represented a better meshing with the finite difference approximations of 
the equilibrium equations than the first element and should, therefore, 
lead to even better agreement with the exact solution. In some respects, 
this was not to be the case. 

Moment equilibrium and force equilibrium for the y direction yield the 
formulas 


r 


m 


m 


lyCcr^Cf^ y) dy]A 


h[a («■ y) dy] + 

xy , 
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Cf (-&» y)dy + 
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y) dy 


» sh. 
2 


0 . 


(B4) 


Using second order Lagrange interpolations for stresses along the edges of 

the material element leads to two formulas that involve the stress 

components (a a a ) at the nine finite difference nodes that are 

X > y > xy 

nearest to the element. Subsequent use of the symmetrical stress-displacement 
relations leads to two equations in the displacements (u,v) . These equations 
are not presented here because of their length, and because no new Insights 
can be derived from them. 

Using only the moment equilibrium equation and the finite difference 
normal stress boundary condition (Similar to the procedure used in the 
analysis of the first material element) it was observed that the stress 
components and agreed well with the exact solution, but was grossly 
over estimated. 

Introducing the force equilibrium equation for the y direction, 
together with the moment equilibrium equation, revealed that a a , a 

were each in substantial disagreement with the exact solution for these 
stress components near the boundary point (a,o) . 

The results of the numerical studies presented in this appendix 
indicate that substantial analytical difficulties are encountered when a 
finite material element is introduced near the point (a,o) to account for a 
lack of S3mimetry in the stress tensor. 

To be analytically rigorous unsymmetrical stress-strain relations should 
be used, and all three equilibrium equations associated with the element 
should be used. 
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FIGURE Bl. Finite material element using linear variation of stresses 
along its edges. 



X - finite difference grid 


0 (x.o) 

y 


points , 


FIGURE B2. Finite material element using quadratic variation of 
stresses along its edges. 
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APPENDIX C 


USER INSTRUCTIONS FOR SEMI-INFINITE PLANE UNDER PARTIAL LOAD. 


This appendix contains information describing the cards that must be 
prepared by the program user and the output information to be expected. 


INPUT INFORMATION 


FIRST CARD (3FX0.0; 415) 


columns 1-10 
11-20 
21-30 
31-35 


36-40 

41-45 

46-50 


Poisson’s ratio, v 

Modules of elasticity, E(psi) 

Uniform pressure, p (Ib/in) 

Column number of the grid line through 
point A, (NDSCN) . NDSCN » 3 in Figure 
C2. 

NROW, number of grid lines parallel to 
the x-axis , NROW « 7 in Figure C2 . 
NCOL, number of grid lines parallel to 
the y-axis. NCOL = 9 in Figure C2. 
NPAR, parameter used to select boundary 
conditions imposed at point A, See 
Boundary Conditions. 


SECOND CARD (6F 10.0) 


Grid line spacings (inches) in the 
x-direction, dx^ are shown in Figure C2. 

THIRD CARD (6F 10.0) 


Grid line spacings (inches) in the 
y-direction, dy^ are shown in Figure C2. 

OUTPUT INFORMATION 


Input data is printed, followed by the stress components a^(SIGX) , 

a (SIGY) , and a (SIGXY) associated with the line x ~ a. Stresses 
y xy ' 

along the loaded boundary are also printed. 


MINIMUM DIMENSIONS FOR ARRAYS 


If program capacity needs to be enlarged the following arrays 
require minimum dimensions as indicated. 
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DIMENSION DX(NCOL-l), DY(NROW-l), ST(2, 18), XST (MA) , 
SIGX(NROW), SIGY(NROW), SIGXY(NROW) , R<NEQ) , 

ID(2, NUMNP), IDIAG(2), ID1(18) , IDD(NROW) , 

Y(NROW) 

NCOL - Number of grid lines parallel to the y-axis. 
NROW - Number of grid lines parallel to the x-ascis, 
NUMNP - Number of nodal points (NR0W*NC0L) , 

. NEQ - 2*NUMNP - NROW - NCOL + 1 

MA - NEQ*(8’^NROW-3)-(2*NROW-l)*(4*NROW-l) 


BOUNDARY CONDITIONS SPECIFIED AT FINITE DIFFERENCE GRID POINTS 
(See Figure C2) . 

f Equilibrium equations 

0 u==v=0 (x and y displacements) 

0 u"0 ,0 =» 0 

xy 

X 0 “0, a =» -p 

xy ’ y 

■ (1) NPAR=1 - a » -p/2, a « 0 

y * xy 

(2) NPAR=2 - -p, finite moment equation. 

(3) NPAR=3 - a^= **p/2, finite moment equation. 

(4) NPAR=4 - Finite force and finite moraf^.nt equations. 

(5) NPAR«5 - a = -p, a »0. 

y ,xy 

(6) NPAR=6 - Oy® ■’P/2, modified moment equation. 

(7) NPAR=7 - Modified finite force and finite moment 

equations . 


A 

□ 

A 


a -a ®0 
y xy 


a -a =0 
X xy 


v=0 , a ®0 
^ xy 


OF POOR QUALITY 
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NOTES : 


(1) H and V in Figure C2 should be large enough to validate the 
assumption that the stresses on the bottom boundary and right 
boundary negligibly small. 

(2) A vertical grid line through point A is required, 

(3) The number of grid lines to the left of point A must be greater 
than two, 

(4) Numbering scheme for nodal points is shown in Figure C2. 
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APPENDIX D 


USER INSTRUCTIONS FOR BIMETALLIC PLATE IN TENSION 


This appendix contains information describing the cards that must 
be prepared by the program user and the output to be expected. 


INPUT INFORMATION 


FIRST CARD (3F10.0, 2l5) 


Columns 1-10 
. 11-20 
21-30 
31-35 

36-40 


Poisson's ratio, v 

Modulus of elasticity, E(psi) 

Uniform pressure, p (Ib/in) 

NROW - Number of grid lines parallel 
to the X axis . 

NCOL - Number of grid lines parallel 
to the y axis 


SECOND CARD (6F10.0) Finite difference spacings for the x 

direction (dx^) (inches), 

THIRD CARD (6F10.0) Finite difference spacings for the 

y direction (dy^) (inches) . 


OUTPUT INF0Ri«IATI0N 


Input data are printed followed by the stress components a^(SIGY) 

and a^y(SIGXY) along the bond line. Stress 

components a^(SIGX) and o^^(SIGXY) are printed for the stress free 

edge. The net normal force and net shear force acting on the bond 
line are also printed. 


MINIMUM DIMENSIONS FOR ARRAYS 

If program capacity needs to be enlarged the following arrays 
require minimum '/imensions as indicated. 

DIMENSION DX(NCOL-l), DY(NROW-l) , ST(2,18), XST(MA) , X(NCOL) , 
Y(NROW), R^EQ), ID(2, NUMNP) , IDIAG(2) , ID1(18) , 
SIGXYL(NROW) , SIGXL(NROW) , SIGXYT(NCOL) , SIGYT(NCOL). 

NCOL - Number of grid lines parallel to the y axis. 
NROW - Number of grid lines parallel to the x axis. 
NUMNP - Number of nodal points (NCOL^NROW) . 

NSQ - Number of equations = 2*(NUMNP-NR0W)-NC0L + 1 
MA - (8*NROW-3)*NEQ-(2*NROW-l)*(4*NROW-l) 
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BOUNDARY CONDITIONS SPECIFIED AT FINITE DIFFERENCE GRID POINTS 
(See Figure D2) 


0 Equilibrium equations. 

X u«v=0 (x and y displacements) 


□ 

a = 

a = 0 

xy 

y 

A 

0 “0 

a =p 


xy 


0 

v=0, 

a =» p/2 

X * 


t v“0 , a =0 
’ xy 


Notes : 

(1) Subroutine DCSPQU is a subroutine from PORT Mathematical 
Subroutine Library and is not provided V7ith, the main program. 
This subroutine is used to integrate the stresses along a 
boundary to determine the net normal force, and net shear 
force on that boundary. If users do not have access to 
PORT subroutines, delete program statements delimited by 
comment statement C /////// in the main program. 

(2) Subroutine DGELB is a subroutine from the IBM Scientific 
Subroutine Package and is included with the main program. 


OF POOR QUALITY 
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APPENDIX E 


USERS INSTRUCTIONS FOR UNIFORM STRAIN OF LAYERED COMPOSITE 

This appendix contains information describing the cards that must be 
prepared by the program user and the output information to be expected. 


INPUT INFORMATION 


FIRST CARD (3F 20.0) 


columns 

1-20 

Poisson's ratio, ('^i2~^23 "^13^ 


21-40 

Shear modulus ^^j^2~^13~^23 ’ 


41-60 

Modulus of elasticity parallel to fiber 
direction, Pascal) 

SECOND CARD 

(2F20.0, 

3110) 

columns 

1-20 

Moduli of elasticity perpendicular to the fiber 
direction ^^22~^33’ • 


21-40 

Applied uniform strain 

< 

41-50 

NROW - Number of grid lines parallel to the y axis. 


51-60 

NCOL - Number of grid lines parallel to the z axis. 


61-70 

INT - Row number for the grid line coincident with 
the interface between the +45 and -45 plies. (INT=5 
in Figure E2) . 

THIRD CARD (6F10.0) 

Grid line spacings (meters) in the y direction, dy. 
are shown in Figure E2. 

FOURTH CARD 

96F10.0) 

Grid line spacings (meters) in the z direction, dz , 
are shown in Figure E2. ^ 


OUTPUT INFORMATION 

Input data are printed, followed by the stress components CT (SIOZ) , 

z 

o^^(SIGZX), o^y(SIGZY), a^(SIGX), a^(SIGY), and cr^y(SIGXY) along the 

interface. Subsequently the stress components and along the 

stress free edge are printed. 


Subscripts u or £. are affixed to a stress component to indicate its 
relation to either the +45 or the -45 ply, respectively. 
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MINIMUM DIMENSIONS FOR ARRAYS 

If program capacity needs to be enlarged the following arrays require 
minimum dimensions as indicated, 

DIMENSION ST(3,27), IDl(£7), ID1AG(3) , 

DY(NCOL-l), DZ(NROW-l), Z(NROW) , 

SGY(NROW), SGXY(NROW), SIGX(NROW) , XST(MA) , 

ID(3, NUMNP), R(NEQ), Y(NCOL) 

NCOL - Number of grid lines parallel to the z axis. 

NROW “ Number of grid lines parallel to the y axis, 
NUMNP- Number of nodal points (NCOL^NROW) 

NEQ - Number of equations (3*NUMNP-2*NR0W - NCOL + 3) 

MA - NEQ*(12*NROW+l)-(3*NROW)*(6*NROW + 1) 


BOUNDARY CONDITIONS SPECIFIED AT FINITE DIFFERENCE GRID POINTS 

(See Figure E2) , Six different programs EDGSTRSl, EDGSTR2 EDGSTRS6 

incorporate the following boundary conditions . 

S U=V=W“o (x, y, and z displacements) 


0 


A 

□ 


A 


X 

o 


A 


U-r:V=W =0 

,y 

a =a =0 =0 

z zx zy 

0 =0 “O =0 

y yx yz 

W=U =V =0 
,z ,z 

Equilibrium equations for the +45 ply. 

Equilibrium equations for the -45 ply. 

Stress continuity along the interface (^ 2 u”^z£J °^zxu~^zx£ 

and a =cr ^ ) 
zyu zy.c 

U=V=W =0 for EDGSTRSl through EDGSTR5 

,y 

Stress continuity between pKes for EDGSTRS6 


0 


0 


W=U. = V =0 for EDGSTRSl 

H 2 j 2 

a =a =0 =0 for EDGSTRS2 through EDGSTR6 

y yx yz “ 

,r Interlaminar stress continuity for EDGSTRSl, EDGSTRS2, and EDGSTRS6 

a =a =0 =0 for EDGSTRS3 

yu yxu yzu 

“ EDGSTRS4 

0,r, +0 p=0„ , +0 /7='0 +0, <5 = 0 for EDGSTRS5 

V. yu y.t yxu yxc yzu yzil 
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NOTES j 


(1) These programs determine the distributions of the stress components 
along the interface between the +45 and -45 plys of a [± 

laminate under uniform axial strain with laminate properties ^ ^22* 


^ 22 “^ 33 * ®12 ■" ^13 “ ^ 3 ’ 


^12 “'^23 


13 


(2) Subroutine DGELB is a subroutine from the IBM Scientific Subroutine 
Package and is provided with the main program. 

(3) Subroutine DCSPQU is a subroutine from PORT Mathematical Subroutine 
Library and is not provided with the main program. This subroutine is 
used to integrate the stresses along a boundary to determine the net 
normal force and net shear forces on that boundary. If users do 

not have access to PORT subroutines, delete program statements 
delimited by comment statement C///// in the main program. 
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